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Abstract 

The  structure  of  receptive  fields  in  the  visual  cortex  is  believed  to  be  shaped  by 
unsupervised  learning  .  A  simple  variant  of  unsupervised  learning  is  the  extraction  of 
principal  components.  In  this  paper,  we  derived  analytically  the  form  of  the  principal 
components  of  natural  images.  This  derivation  relies  on  results  about  the  covariance 
matrix  of  natural  images  (Field,  1987).  Our  resluts  predict  both  the  shapes  and  the 
phases  of  the  receptive  fields.  We  also  compared  our  resluts  to  numerical  simulation 
results  (Hancock  et  al.,  1992).  Finally  the  biological  relevance  of  our  results  is  discussed. 


1  Introduction 

It  is  generally  believed  that  the  shape  of  the  receptive  fields  in  the  visual  cortex  is  determined 
by  some  form  of  unsupervised  learning.  Hebb’s  postulaie  (Hebb,  1949)  is  the  cornerstone 
of  most  unsupervised  learning  models  in  neural  networks.  The  naive  Hebbian  rule,  how¬ 
ever,  is  unstable.  This  problem  is  often  overcome  by  adding  constraints  to  this  rule.  One 
such  constraint  is  adding  a  decay  term  to  this  rule.  It  has  been  proven  that  a  neuron  wulh 
Hebbian  learning  rule  plus  a  proper  decay  term  can  perform  a  principal  component  extrac- 
"tion  (Oja,  1982).  Furthermore,  a  neural  network  with  proper  lateral  inhibition  can  perform 
the  extraction  of  several  principal  components  simultaneously  (Sanger,  1989).  The  compu¬ 
tational  importance  of  principal  components  is  that  they  are  the  optimal  linear  projections 
for  minimizing  the  mean  squared  reconstruction  error  (Fukunaga,  1990). 
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Since  the  principal  components  of  a  set  of  inputs  depend  only  on  their  covariance  ma¬ 
trix,  it  is  reasonable  that  given  this  matrix,  they  can  be  calculated  analytically.  A  known 
example  is  Linsker’s  multilayered  network  (Linsker,  1986;  Linsker,  1988),  which  is  trained 
with  random  noise.  The  covariance  matrix  results  from  the  existence  of  the  Arbor  functions 
between  the  random  input  layer  and  the  first  layer.  This  induced  correlation  was  used  by 
several  investigators  to  give  an  analytic  explanation  of  Linsker’s  result  (Kammen  and  Yuille, 

1988;  Yuille  et  ah,  1989;  Tang,  1989;  MacKay  and  Miller,  1990a;  MacKay  and  Miikr,  1990h) 
1 

We  believe  that  modeling  the  environment  with  natural  scenes  is  more  reasonable  than 
modeling  it  with  random  noise  or  bars  or  edges,  etc.,  because  the  environment  to  which  the 
visual  cortex  is  exposed,  through  the  visual  pathway,  is  composed  of  natural  scenes.  When 
the  input  environment  is  not  random,  the  properties  of  the  covariance  matrix  depend  on 
the  nature  of  the  environment,  and  on  the  preprocessing  performed  by  the  visual  pathw’ay. 
However,  in  this  paper,  we  have  simplified  the  problem  by  neglecting  the  effects  of  this 
preprocessing,  as  a  first  step  towards  a  more  complete  mathematical  understanding  of  the 
receptive  field  structure  of  neurons  in  the  visual  cortex.  Another  .aspect  of  this  simplification 
is  that  we  can  use  our  results  to  explain  the  simulation  results  (Hancock  et  ah,  1992),  in 
which  the  authors  used  a  network  of  tne  type  proposed  by  Sanger  (Sanger,  1989)  to  extract 
principal  components  from  natural  scenes.  In  their  simulations,  no  preprocessing  was  applied 
to  the  images,  apart  from  multiplying  the  inputs  by  a  smooth  Gaussian  function,  in  order 
to  eliminate  edge  effects. 

The  nature  of  the  covariance  matrix  ol  natural  images  was  investigated  by  Field  (Field, 
1987).  He  found  that  the  spectrum  of  covariance  matrix  is  proportional  to  the  inverse 
of  the  square  of  the  frequency.  In  section  2,  we  will  derive  analytically  the  form  of  the 
principal  components  of  natural  images  under  this  assumption..  Rather  t  han  taking  a  smooth 
Gaussian  window  as  Hancock  et.  al.  did,  we  aissume  a  circular  hard  boundary  to  the  receptive 
fields.  We  find  that  the  solutions  are  the  Fourier- Bessel  functions  (Jackson,  1975,  etc.).  We 
will  show  in  section  3,  that  under  the  assumption  that  the  covariance  matrix  spectrum  has 
a  small  non-rotationally  symmetric  correction,  the  solutions  have  a  definite  phase. 

When  simulation  results  (Hancock  et  al.,  1992)  are  compared  to  our  results,  we  find  a 
good  agreement  for  the  higher  eigenvalue  solutions,  with  some  deviations  for  low'er  eigenvalue 
solutions.  In  the  discussion  we  explain  how  these  deviations  may  come  about. 

It  is  clear  that  the  receptive  fields  obtained  here  are  not  identical  to  the  receptive  field 
structures  obtained  experimentally  in  the  visual  cortex.  Does  this  imply  that  neurons  in  the 
first  layers  of  the  visual  cortex  are  not  principal  component  analyzers?  We  will  address  this 
issue  in  the  discussion. 

*Their  solutions  are  to  an  equation  with  some  additional  terms  to  the  principal  component  equation,  but 
it  is  similar  to  a  principal  component  equation. 
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2  The  Rotationally  Symmetric  Solution 

The  principal  components  are  the  eigen-functions  of  the  covariance  matrix.  Therefore  the 
equation  we  try  to  solve  is  the  eigenvalue  problem,  i.e.,  the  eigen-equation,  which  has  the 
form 

X!  =  Au;,  ( 1 ) 

1 

where  w,  are  eigen-vectors,  A  is  the  eigenvalue,  and  C,j  is  the  covariance  matrix  which  is 
defined  as  C.j  =  £[{1^  —  £![/,])(/j  —  £^[/_,])]  for  input  pattern  {/,}•  Since  we  are  dealing  with 
two  dimensional  space,  the  index  i  really  denotes  a  point  in  the  tvjo  dimensional  space,  so  it  is 
more  convenient  to  rewrite  the  covariance  matrix  in  the  form  C(r,,  rt).  Due  to  translational 
invariance,  C'(r,,r')  =  C(r,  —  rf).  In  the  continuous  limit,  the  summation  will  become  an 
integral  o'''er  r\  th'’«  the  eigen-equation  becomes 

j  C{r  —  r')iy(r')d^r'  =  \w{r).  (2) 


in  which  t(;(r)  is  the  continuous  limit  of  the  eigen- vectors  u;,. 

The  Fourier  transform  (spectrum)  of  the  covariance  matrix  has  the  form,  C'(k)  --  cjV} 
where  c  is  a  constant  (Field,  1987).  Hereafter  w'e  wrill  set  c  —  1  for  convenience.  Thus  C(r) 
satisfies  V^C'(r)  =  — ^(r)  which  can  be  readily  proven  by  taking  Fourier  transformation  on 
both  side  of  this  equation.  Therefore  by  operating  both  sides  of  the  eigen-equation  2  with 
we  get 

V^ro(r)  =  Tti^(r).  (3) 

A 

In  order  to  find  a  closed  form  solution  to  this  equation,  w'e  assume  that  w[r)  is  non-zero 
only  within  a  circle  of  radius  a.  The  justification  for  this  is  that  the  receptive  field  of  a 
neuron  is  always  of  finite  size.^  The  equation  can  be  expressed  in  the  polar  coordinate,  and 
the  solutions  are  (Jackson,  1975,  etc.) 


V 


co5(m^) 

sin{m9) 


for  r  <  a 
for  r  >  a 


(4) 


in  which  m  =  0,  1,  2,...,  Jmi^)  is  the  standard  Bessel  functions,  A„„  is  the  rth  root  of 
equation  Jm{ajy/X)  =  0  ,  r  and  6  are  the  polar  coordinate  of  r.  Because  the  complete 
solutions  of  the  original  eigenvalue  problem  must  be  a  complete  and  orthonormal  set  of 
functions,  the  complete  set  of  solutions  to  the  eigenvalue  problem  is  ^ 

Gaussian  mask  used  in  the  numerical  simulation  (Hancock  et  al.,  1992)  is  a  relaxed  version  of  the 
boundary  condition  we  have  here.  Solving  the  exact  equation  with  a  Gaussian  mask  becomes  mathematically 
more  complicated.  We  should  expect  that  for  a  Gaussian  mask  the  solutions  near  the  center  should  be  almost 
identical  to  the  hard  boundary  solutions,  but  they  would  approach  zero  gradually  near  the  boundary. 

^Of  course  the  trivial  solution  that  wjl,j(r)  equals  zero  must  be  deleted  from  this  set. 
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w  N  (  •/m(-7f-)co5(m0  4-  for  r  <  a 

wL  ( r)  =  s 

\  0  for  r  >  a 

2  /  ^  f  Jm{-i^)sin{Tn6  +  4>r^,)  for  r  <  a 

^  ^  I  0  for  r  >  a 


(5) 


where  is  a  set  of  undetermined  phases.  These  two  eigen-functions  have  the  same  eigen¬ 
value  Xmi,  i-e.,  they  are  degenerate. 

If  we  order  the  solutions  by  the  magnitudes  of  the  correspondent  eigenvalues  A^,,  the 
first  ten  solutions,  with  4>mi  ~  0  and  a  =  1,  are  drawn  in  figure  1. 

3  Retrieving  the  Phase 

The  solutions  above  and  w^,(r)  not  only  have  undetermined  phases,  but  also  are 

degenerate.  This  contradicts  the  results  of  the  simulations  (Hancock  et  al.,  1992)  in  which 
the  phases  seem  to  always  take  the  value  zero,  and  the  solution  has  a  different  eigenvalue 
from  the  solution.  These  results  can  be  retrieved  if  we  assume  that  the  covariance  matrix 
has  a  non-rotationally  symmetric  perturbation  term.  This  assumption  is  not  arbitrary  since 
an  inspection  of  Fields  results(Field,  1987,  figure  7)  reveals  that  this  is  indeed  the  case. 
Hereafter  we  assume  this  perturbation  term  has,  in  k  space,  the  form 


C'(k)  =  U(k)T{h). 


(6) 


In  order  to  calculate  this  perturbation,  the  representation  of  this  perturbation  in  the  two 
degenerate  eigen- functions  and  has  to  be  calculated.  It  is  easier  to  perform 
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this  in  k  space  in  which  the  eigen-functions  and  n;^,(r)  are  replaced  by  their  Fourier 

transforms, 


~  -r  (f>Tni') 

tafntfk)  "h  ^mt) 


in  which 


fmt{k)  =  ■nj"'  f  Jfn{-  Z  )Jm{kr)rdT 

Jo  VXm. 

where  —  —  1.  If  we  denote 


(7) 

(8) 


r(^k)  =  E 

i 


(9) 


which  is  the  Fourier  expansion  of  T{6^^).  The  representation  of  the  perturbation  matr  ix  '.vi*  h 
respect  to  the  two  degenerate  eigen-functions  has  the  form 


(cos(6)  sin{6)  \ 
sm{6)  -cos{6)  j 


(10) 


in  which  ^  =  2(f)mi-\-2ma2m  &nd  g^i  ~  ^t2m  J  U{k)\fmt{k)\'^kdk.  Since  the  two  eigen-functions 
are  degenerate,  any  linear  combination  of  these  two  eigen-functions  is  an  eigen-function  of  C. 
Therefore,  all  we  have  to  do  is  to  find  a  linear  combination  of  them  which  diagonalizes  the 
perturbation  matrix,  i.e.,  to  find  the  eigenvalues  and  eigen-vectors  of  the  matrix  in  equation 
10,  which  are 


/  cos{8l2)  \ 
sm(6/2)  j 

(11) 

/  — sm(^/2)  \ 
cos{6/2)  j 

with  eigenvalues  gmi  and  —gmx  respectively.  Furthermore,  \iU{k)  —  ejk^  then  the  complete 
expression  for  the  correction  to  the  eigenvalue  takes  the  form  gmi  —  f2m/2. 

Thus  the  eigen-functions  and  eigenvalues  after  the  perturbation  can  be  readily  written 
out  as 


,(k)  :=  J„(  ^  )cos(m(0 

V  '^mi 

-  OL2m)) 

T 

(k)  -  J^{—^--)sin{m[6  - 

0^2m)) 

(12) 
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with  eigenvalues  A  +  ,  =  A^,  +  ,  and  A„,  =  A^,  -  Qm,  ,  respectively.  So  the  degeneracy  is 

broken.  This  is  in  agreement  with  Hancock’s  simulations.  These  solutions  have  an  important 
feature,  i.e.,  their  phases  are  determined  by  the  properties  of  the  real  w'orld  covariance 
matrix.  If  the  covariance  matrix  has  a  definite  symmetry  with  an  inclination  angle  a,  then 
the  solutions  would  also  have  the  same  symmetry  angle.  Because  in  this  case  a2m  =  a  for 
all  m.  The  spectrums  of  the  covariance  matrix,  shown  in  hgure  7  of  Field’s  paper,  indeed 
indicates  a  symmetry  axis  along  a  =  0.  Thus  equation  12  predicts  the  zero  phase  result 
found  in  Hancock’s  simulation.  When  Hancock  used  images  which  were  tilted  by  45  degrees 
before  being  scanned,  the  preferred  axis  of  the  receptive  fields  was  found  to  be  45  degrees. 
Again  this  is  predicted  by  equation  12.  because  the  symmetry  axis  of  the  covariance  matrix 
spectrum  also  gets  rotated  by  45  degrees  due  to  the  rotated  images,  i.e.,  a  =  45°,  and  thus 
the  solutions  also  get  rotated  by  45  degrees. 

4  Discussion 

We  have  calculated  the  forms  of  the  principal  components  of  natural  images  based  on  the 
result  about  the  covariance  matrix  (Field,  1987),  and  have  shown  that  a  non-rotaticnal!y 
symmetric  perturbation  can  break  the  degeneracy  and  give  a  definite  phase  which  only 
depends  on  the  properties  of  the  real  world  covariance  matrix.  These  results  for  a  large  part 
agree  with  the  numerical  simulation( Hancock  et  al.,  1992). 

A  discrepancy  can  be  noticed  between  our  results  and  Hancock’s  simulation  results  for 
eigen-functions  corresponding  to  smaller  eigenvalues.  This  comes  from  the  mixing  of  solu¬ 
tions,  due  to  a  perturbation  added  to  the  covariance  matrix(Merzbacher,  1970).  This  mixing 
would  be  more  noticeable  for  lower  eigenvalue  solutions  since  for  smaller  A„,  and  i.e., 
for  large  m  and  m'.  Ami  ^  •^m'i,  the  solutions  are  nearly  degenerate.  Thus  any  small  per¬ 
turbation  will  yield  highly  mixed  linear  combinations  of  different  eigen-functions,  the  same 
linear  combination  phenomena  we  have  discussed  in  section  3.  The  perturbation  can  come 
from  the  difference  between  the  Gaussian  mask  used  in  the  simulation,  and  the  hard  bound¬ 
ary  we  chose  in  the  analytic  analysis.  It  can  also  come  from  all  sorts  of  correction  to  the 
spectrum  of  the  covariance  matrix  including  the  non-rotationally  symmetric  correction.  An¬ 
other  correction  can  be  a  high  frequency  cutoff  due  to  the  finite  sampling  frequency  of  the 
images. 

The  neurobiological  relevance  of  the  type  of  technique  used  in  this  paper  is  that  we 
can  deduce  for  different  learning  rules  what  kinds  of  receptive  fields  they  should  produce. 
Given  these  receptive  fields,  we  can  compare  them  to  the  real  biological  receptive  fields. 
This  comparison  can  be  used  to  assess  whether  the  biological  hardware  really  implements  or 
approximates  a  theoretically  proposed  learning  rule. 

The  most  obvious  conclusion  which  stands  out  when  we  observe  the  results  in  figure  1 ,  is 
that  these  receptive  fields  have  little  resemblance  to  receptive  fields  reported  in  the  biological 
literature  (Hubei  and  Wiesel,  1959,  etc.).  Does  this  imply  that  biological  neurons  are  not 
principal  component  analizers?  When  addressing  this  question  we  have  to  keep  in  mind  that 
the  natural  images  projected  on  the  retina,  undergo  preprocessing  in  the  retina  and  LGN, 


before  they  reach  the  visual  cortex.  Similar  preprocessing  should  therefore  be  applied  to 
natural  images  in  simulations  and  analytic  studies,  before  a  sensible  answer  can  be  given. 
What  can  also  be  observed  is  that  the  receptive  fields  of  the  neurons  in  visual  cortex  seem  to 
form  a  non-orthogonal  set  such  as  Gabor  function  family  (Daugman,  1985,  etc.)  as  opposed 
to  the  orthogonal  and  complete  function  set  we  got  here.  This  may  result  from  neurons 
being  noisy  entities,  a  property  of  neurons  not  being  taken  into  account  when  formulating  the 
minimal  mean  square  reconstruction  error  criterion,  which  results  in  the  principal  component 
eigen-equation.  All  these  issues  are  currently  under  study. 
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